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^ABSTRACT: In the field of structural reliability, the Monte-Carlo estimator is considered as the reference 
^probability estimator. However, it is still untractable for real engineering cases since it requires a high number 
5_iOf runs of the model. In order to reduce the number of computer experiments, many other approaches known 
PLas reliability methods have been proposed. A certain approach consists in replacing the original experiment by 
"^a surrogate which is much faster to evaluate. Nevertheless, it is often difficult (or even impossible) to quantify 
QQthe error made by this substitution. In this paper an alternative approach is developed. It takes advantage of 
^the kriging meta-modeling and importance sampling techniques. The proposed alternative estimator is finally 
I applied to a finite element based structural reliability analysis. 



^ 1 INTRODUCTION 

I ^Reliabilitv analysis consists in the assessment of the 
level of safety of a system. Given a probabilistic 
^ model (a random vector X with probability density 
s^ofunction (PDF) /) and a performance model (a func- 
I^tion g), it makes use of mathematical techniques in 
"border to estimate the system's level of safety in the 
^form of a failure probability. A basic approach, which 



limit theorem, this estimator is asymptotically unbi- 
ased and normally distributed with variance: 
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In practice, this variance is compared to the unbiased 
estimate of the failure probability in order to decide 
whether it is accurate enough or not. The coejficient 



^ makes reference, is the Monte-Carlo simulation tech- of variation is defined as 
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^nique that resorts to numerical simulation of the 
^performance model through the probabilistic model. 
^Failure is usually defined as the event F ^{g {X) < 0}, 
•^so that the failure probability is defined as follows: 

Cipf^¥{F)=P{{g{X)<0})=[ f{x)dx (1) 



Introducing the failure indicator function lg<o being 
equal to one if g (x) < and zero otherwise, the failure 
probability turns out to be the mathematical expecta- 
tion of this indicator function with respect to the joint 
probability density function / of the random vector 
X. The Monte-Carlo estimator is then derived from 
this convenient definition. It reads: 
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where ... e N* is a set of samples 

of the random vector X. According to the central 



N, this coefficient dramatically increases as soon as 
the failure event is too rare (pj ^ 0) and proves that 
the Monte-Carlo estimation technique intractable for 
real world engineering problems for which the perfor- 
mance function involves the output of an expensive- 
to-evaluate black box function - e.g. a finite element 
code. Note that this remark is also true for too fre- 
quent events (pj 1) as the coefficient of variation of 
I- PfMC exhibits the same property. 

In order to reduce the number of simulation runs, 
a large set of other approaches known as reliability 
methods have been proposed. They might be classi- 
fied as follows. 

A first approach consists in replacing the orig- 
inal experiment by a surrogate which is much 
faster to evaluate. Among such approaches, there 
are the well-known first and second order re- 
liability methods (e.g. [Ditlevsen and Madsen 



'1996; 'Lem aire[ 2009[), q uadratic response surfaces 
(Bucher and Bourgund[ [T990) and the more re- 
cent meta-models such as support vector machines 



( |Hurtado[ |2004[ [Deheeger and Lemaire] |2007|), neu- 
ral netw orks ([Pap adrakakis and Lagaros', '2002) and 



kriging (|Kaymaz 



.2005^ ,Bichon et al^ |2008). Never- 



theless, it is often difficult or even impossible to quan- 
tify the error made by such a substitution. 

The other approaches are the so-called variance re- 
duction techniques. In essence, these techniques aims 
at favoring the Monte-Carlo simulation of the fail- 
ure event F in order to reduce the estimation vari- 
ance. These approaches are more robust because they 
do not rely on any assumption regarding the func- 
tional relationship g, though they are still too compu- 
tationally demanding to be implemented for industrial 
cases. For an extended review of these techniques. 



the reader is referred to the book by Rubinstein and 



Kroese(2008). 



In this paper an hybrid approach is developed. It 
is based on both margin meta-models (defined here- 
after) and the importance sampling technique. It is 
then applied to an academic structural reliability prob- 
lem involving a linear finite-element model and a two- 
dimensional random field. 

2 ADAPTIVE PROBABILISTIC 

CLASSIFICATION USING MARGIN 
META-MODELS 

A meta-model means to a model what the model it- 
self means to the real-world. Loosely speaking, it is 
the model of the model. As opposed to the model, its 
construction does not rely on any physical assump- 
tion about the phenomenon of interest but rather on 
statistical considerations about the coherence of some 
scattered observations that result from a set of experi- 
ments. This set is usually referred to as a design of ex- 
periments (DOE): X = {xi, ... , Xm}. It should be care- 
fully selected in order to retrieve the largest amount of 
statistical information about the underlying functional 
relationship over the input space Vx. Here, we attempt 
to build a meta-model for the failure indicator func- 
tion lg<o. In the statistical learning theory (jVapnik 



1995 ) this is referred to as a classification problem. 

Hereafter, we define a margin meta-model as a 
meta-model that is able to give a probabilistic predic- 
tion of the response quantity of interest whose spread 
{i.e. variance) depends on the lack of information 
brought by the DOE. It is thus reducible by bring- 
ing more observations into the DOE. In other words, 
this is an epistemic (reducible) source of uncertainty. 
To the authors' knowledge, there exist only two fam- 
ilies of such margin meta-models: the probabilistic 
support vector machines P-SVM by |Platt| ( |1999p and 
Gaussian-Process- (or kriging-) based classification 



( Santner et al. 2003 1. The present paper makes use of 
the kriging meta-model, but the overall concept could 
easily be extended to P-SVM. The theoretical aspects 
of the kriging prediction are briefly introduced in the 
following subsection before it is applied to the classi- 
fication problem of interest. 



2.1 Gaussian-process based prediction 

The Gaussian-Process based prediction (also known 
as kriging) theory is detailed in the book by Santner 
et al. ( 2003| ). In essence, kriging assumes that the per 



formance function g is a sample path of an underlying 
Gaussian stochastic process G that would read as fol- 
lows: 



G{x) = f{xfp + Zix) 
where: 



(4) 



• / (cc)^ f3 denotes the mean of the GP which cor- 
responds to a classical linear regression model 
on a given functional basis {/j, i = ... ,p} e 
C2 iVx,R) ; 

• Z (x) denotes the stochastic part of the GP which 
is modelled as a zero mean, constant variance o-^, 
stationary Gaussian process with a given sym- 
metric positive definite autocorrelation model. It 
is fully defined by its autocovariance function 
which reads ((a;, x') eV^ x Vx). 

Cgg{x,x') ^aQR{\x-x\,e.) (5) 

where £ is a vector of parameters defining R. 

The most widely used class of autocorrelation func- 
tions is the anisotropic squared exponential model: 



R(^\x — x'\ ,£) = exp [ ^ - 
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The be st linear unbiased est i mation of G at p oint x 
is shown ( [Santner et alj |2003| |Severini[ |20051 Chap. 
8) to be the following Gaussian random variate: 



G{x) = G{x) \{g{xi), ... ,g{xm)} 
with moments: 

//g [x) = f {xf ^ + r {xf R-Uy-F$ 



(7) 
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where we have introduced r, R and F such that: 



ri{x) ^ R{\x - Xi\, £) , i = I, . . 

J — ^ I % J I ' ^ ) ' ^ — ^ ' 



, m 

. . ,m, j = 1, 
1, . . . ,m 



(10) 

,m (11) 
(12) 



At this stage it can easily be proven that /ig (xi) = 
g (xi) and (Tg (aji) for i = 1, ... , m, thus meaning the 
kriging surrogate is an exact interpolator. 



Given a choice for the regression and correlation 
models, the optimal set of parameters (3*, £* and 
(T^* can then be inferred using the maximum likeli- 
hood principle applied to the unique sparse obser- 
vation of the GP sample path grouped in the vector 
y = {g (xi) , i = 1, . . . , m) . This inference problem turns 
into an optimization problem that can be solved an- 
alytically for both /3* and a^* assuming £* is known. 
Thus, the problem is solved in two steps: the max- 
imum likelihood estimation of £* is first solved by a 
global optimization algorithm which in turns gives the 
optimal values of /3* and a^f. 

2.2 Probabilistic classification function 

A surrogate-based reliability analysis simply consists 
in replacing the performance function g by its meta- 
model fiQ. Again, this meta-model may be a first- 
or second-order Taylor expansion of the limit-state 
surface g = at a so-called design point ([Ditlevsen 



[and Madsen} [T996| FORM/SQRM), a polyno mial re 
sponse surface ( [Bucher and Bourgund||1990|), a neu- 
ral networks based prediction ( [Papadrakakis and La- 
garos, 2002), or a kriging based prediction ( [Bichon 



et al. 2008 [ )! This surrogate may not be fully accurate 



and it is difficult or even impossible to quantify the 
error made by substitution on the final failure proba- 
bility of interest. 



In this paper as in the work by Picheny (2009), 
it is proposed to use the complete probabilistic pre- 
diction provided by the kriging meta-model instead 



of the sole mean prediction (e.g. as in Bichon et al. 



2008 ). Indeed, since the probabilistic distribution of 



the prediction is fully characterized, the probability 
that the prediction is negative may be expressed in 
closed-form and reads as follows: 



Gix) < 
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Note that this latter quantity is not the sought failure 
probability pj, this is simply the probability that the 
prediction G at some deterministic vector x is nega- 
tive. 

Picheny (|2009') proposes then to use this proba- 



bilistic classification function as a surrogate for the 
real failure indicator function, and uses crude Monte- 
Carlo simulation to estimate the failure probability. A 
different use is proposed in the next section. 

Figure [T] illustrates the concepts introduced in this 
section on a basic structural reliability example from 
Der Kiureghian and Dakessian' (1998). This exam- 
ple involves two independent standard Gaussian ran- 
dom variates Xi and X2, and the performance function 
reads: 



g{xi, X2) ^ b - X2 - k{xi - eY 



(14) 



where & = 5, k = 0.5 and e = 0.1. In subfigure |l(a) 
the limit-state function g{x) = is represented by the 



black dash-dot line. The red minusses (g < 0) and 
blue plusses (g > 0) represent the initial DOE from 
which the kriging meta-model is built. The mean 
prediction's limit-state ^g(a;) is represented by 
the dashed black line. It can be seen that the meta- 
model is not fully accurate since the green triangle x^ 
(among others) is misclassified. Indeed, x^ is safe ac- 
cording to the real performance function g, but it fails 
according to the mean prediction of the meta-model 
/^g. The probabilistic classification function makes a 
smoother decision possible: x^ fails with a 60% prob- 
ability w.r.t. the epistemic uncertainty in the random 
prediction G{x^) ^ 7V(^g(a;°), (Tg(a;°)). Note also that 



the red and blue points in the DOE fails with prob- 
abilities 100% and 0% (safe) respectively due to the 
interpolating property of the kriging metamodel. Sub- 



figure 1(a) is the one-dimensional illustration of the 
three classification strategies for the vector x^. The 
deterministic decision function is an heaviside func- 
tion centered in zero, and the probabilistic classifica- 
tion is a smoother Gaussian cumulative density func- 
tion. 

2.3 Refinement of the probabilistic classification 
function 

In this subsection, a strategy is proposed in order to 
refine the probabilistic classification function so that 
it tends towards the real indicator function lg<o. 

First, let the margin of uncertainty M be defined as 
follows: 

M= jx : -fccrg(a;) < G(a;) < +fc(Tg(a;)| (15) 

where k might be chosen as fc = <f>~^(97.5%) = 1.96 
meaning a 95% confidence interval onto the predic- 
tion of the limit-state surface is chosen. Such a 95% 
confidence margin is illustrated in Subfigure |l(a) as 
the area bounded below by the blue line (2.5% confi- 
dence level) and above by the red line (97.5% confi- 
dence level). The points that are located in this margin 
have an uncertain sign, the others being either failed 
or safe with a confidence level greater than 97.5%. 

We also define the probability that a point x eV^ 
belongs to this margin of uncertainty. Due to the 
Gaussian nature of the prediction, this probability 
may also be expressed in closed-form and reads as 
follows: 
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ag(x) ) 

-ka-Q (x) - fiQ (x) 
aQ{x) 



(16) 



Then, finding the point that maximizes this quan- 
tity on the support of the PDF of X will finally bring 
the best improvement point in the DOE. Starting with 
this statement, many authors in the kriging literature 
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(a) Limit-states and probability contours. (b) Classification given x = a;°. 

Figure 1: Comparison of the three classification strategies on a two-dimensional example from |Der Kiureghian and Dakessian| ( [T998| ). 



decide to use global optimization algorithms in or- 
der to find the best improvement point. For instance, 
Bichon et al.| (|2008) use a different criterion named 
the expected feasibility function, and |Lee and Jung 
( 2008| ) use the constraint boundary sampling crite- 
rion. Note also that an equivalent concept is used by 
Hurtado' ('2OT); |Deheeger and Lemaire] (|2007|); pe^ 



heeger (2008) ; Bourinet et al.|(|2010| ) for SV M. 

In this paper, as in pubourg et al.| ( |201 1| ), a slightly 
different strategy is proposed in order to add several 
points in the DOE. The proposed criterion P [a; e M] is 
multiplied by a weighting density function w so that 

C{x) ^V{x eM) w{x) (17) 

can itself be regarded as a PDF up to an unknown but 
finite normalizing constant. The weighting density w 
can either be chosen as the original PDF of X, or, as 
it is proposed here, the uniform PDF on a sufficiently 
large confidence region of the original PDF. Such a 
confidence region might be difficult to define for any 
given PDF, but as it is usually done in structural relia- 
bility ( [Ditlevsen and Madsen[|1996| ), the original ran- 
dom vector X can be transformed into a probabilisti- 
cally equivalent standard Gaussian random vector U 
for which the confidence region is simply an hyper- 
sphere with radius Pq. The reader is referred to|Lebrun" 



and Dutfoy (2009) for a recent discussion on such 
mappings U = T{X). In that given space, Pq can be 
easily selected as e.g. /^o = 8 which corresponds to the 



maximal generalized reliability index (Ditlevsen and 



Madsen, 1996) that can be justified numerically, and 
the sought uniform PDF is simply defined in terms of 
the following indicator function: 

Markov-chain Monte-Carlo simulation techniques 
{e.g. the slice sampling technique proposed by Neal, 



2003 [ ) might be used in order to generate N (say 
N = 10^) samples from the pseudo-PDF C. These sam- 
ples are expected to be highly concentrated around 
the maxima of the criterion P(a; e M), and thus in the 
vicinity of the predicted limit-state hq{x) = where 
the sign of G is the most uncertain. This large candi- 
date population can then be reduced to a smaller one 
that condensate its statistical properties by means of 



1967). 



a if-means clustering algorithm ( [MacQueen] 
The K (K being given) cluster centers uniformly span 
the margin M and may be added to the DOE in order 
to enrich the prediction of the performance function 
in the vicinity of the limit- state and thus reduce the 
margin of uncertainty. 



META-MODEL-BASED IMPORTANCE 
SAMPLING 



Picheny ( |2009 ) proposes to use the probabilistic clas- 
sification function as a surrogate for the real indica- 
tor function, so that the failure probability is rewritten 
from its definition in Eq. ([!]) as follows: 



Pfe 



G{x)<0 f{x)dx = Ef P G{X)<0 (19) 



It is argued here that this latter quantity does not 
equal the failure probability of interest because it 
sums the aleatory uncertainty in the random vector X 
and the epistemic uncertainty in the prediction G. This 
is the reason why pf^ will hereafter be referred to as 
the augmented failure probability. As a matter of fact, 
even if the epistemic uncertainty in the prediction can 
be reduced {e.g. by enriching the DOE as proposed in 
section [23]), it is impossible to quantify the contribu- 



tion of each source of uncertainty a posteriori. 



This remark motivates the approach introduced 
in this section where the probabilistic classification 
function is used in conjunction with the importance 
sampling technique in order to build a new estimator 
of the failure probability. 

3 . 1 Importance sampling 



According to [Rubinstein and Kroese (2008), impor- 
tance sampling (IS) is the most efficient variance re- 
duction technique. This technique consists in comput- 
ing the mathematical expectation of the failure indica- 
tor function according to a biased PDF which favors 
the failure event of interest. This PDF is called the 
instrumental density. 

Given an instrumental density h, such that h dom- 
inates lj<Q f, the definition of the failure probability 
of Eq. ^ may be rewritten as follows: 
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h{X) 



(20) 



where the subscript h on the expectation operator is 
added to recall that X is therefore distributed accord- 
ing to h. Note that the domination requirement of h 
over lp<o / simply means that: 

SO that the so-called likelihood ratio £{x) = f{x)/h{x) 
is finite for any given x eVx. 

The latter definition of the failure probability easily 
leads to the establishment of the importance sampling 
estimator which reads as follows: 
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where |a;(^), . . . e N* is a set of samples from 

h. According to the central limit theorem, this estima- 
tion is unbiased and its quality may be measured by 
means of its variance of estimation which reads: 
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Rubinstein and Kroesej ( |2008| ) show that this vari- 
ance is zero (optimality of the IS estimator) when the 
instrumental PDF is chosen as: 



X ^ '^g<oix)f{x) ^ lg<o{x)f{x) 

Jlg<0{x)fix)dX Pf 



(24) 



However this instrumental PDF is not implementable 
in practice because it involves the sought failure prob- 
ability in its denominator. There exists infinitely many 
PDF h that allows to significantly reduce the variance 
of estimation though. 



3.2 A meta-model-based approximation of the 
optimal instrumental PDF 

Different strategies have been proposed in order to 
build quasi-optimal instrumental PDF suited for spe- 
cific estimation problems. For instance, Melchers] 
( |1989| ) uses a standard normal PDF centered onto the 
most probable failure point (MPFP) in the space of 
the independent standard Gaussian random variables 
U = T{X) in order to estimate a failure probability. 
Although this approach may lose accuracy as soo n as 
the MPFP is not unique. Cannamela et al. ( 2008[ ) use 
a kriging prediction of the performance function g in 
order to build an instrumental PDF suited for the es- 
timation of extreme quantiles of the random variate 
G^9{X). 

Here, it is proposed to use the probabilistic classi- 
fication function in Eq. ( [T3] ) as a surrogate for the real 
indicator function in the optimal instrumental PDF in 



Eq. (24). The proposed quasi-optimal PDF thus reads 
as follows: 



h*{x) = 
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where p/e is the augmented failure probability which 
has been already defined in Eq. (19). This quasi- 
optimal instrumental PDF is compared to its impracti- 
cal optimal counterpart in Figure |2] using the example 
of Section [121 

3.3 The meta-model-based importance sampling 
estimator 

Choosing the proposed quasi-optimal instrumental 
PDF in Eq. ( [25] ) in the importance sampling defini- 
tion of the failure probability in Eq. ( [20] ) leads to the 
following new definition: 



Pf= I lg<0ix)M^l?{x)dx 

h*[x) 



■Pfe 



Pfe' 



G{x) < 



h* (x) dx 



(26) 



(27) 



(28) 



where we have introduced: 



acorr = Er 
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(29) 



This means that the failure probability is now defined 
as the product between the augmented failure proba- 
bility and a correction factor aeon- This correction 
factor is defined as the expected ratio between the real 
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(a) The optimal instrumental PDF. 



(b) A quasi-optimal PDF. 



Figure 2: Comparison of the instrumental PDF on the two-dimensional example from Der Kiureghian and Dakessian ( 1998 1. 



indicator function lc,<o and the probabilistic classifi- 
cation function P[G(») < 0]. Thus, if the kriging pre- 
diction is fully accurate, the correction factor is equal 
to unity and the failure probability is identical to the 
augmented failure probability (optimality of the pro- 
posed estimator). On the other hand, in the more gen- 
eral case where the kriging prediction is not fully ac- 
curate, the correction factor modifies the augmented 
failure probability accounting for the epistemic uncer- 
tainty in the prediction. 

The two terms of the latter definition of the failure 
probability may now be estimated using Monte-Carlo 
simulation: 
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(30) 
(31) 



where the first iV^- sample is generated from the orig- 
inal PDF /, and the second A^con-sample is generated 
from the quasi-optimal instrumental PDF h*. Accord- 
ing to the central limit theorem, these two estimates 
are unbiased and normally distributed. Their respec- 
tive variance of estimation denoted by erf and Ccorr 
not given here but they might be easily derived. 

To generate samples from h*, it is proposed to 
use a Markov chain Monte-Carlo simulation tech- 
nique which is applicable to a broad class of improper 
PDF for which the normalizing constant is not known 
and thus for the instrumental PDF of interest h*{x) oc 
V[G{x) < 0] f{x). The work presented here makes use 
of the slice sampling technique ( Neal[ 2003 1. 

Finally, the final estimator of the failure probability 
simply reads as follows: 



The calculation of the coefficient of variation of the fi- 
nal estimator pymetais will be detailed in a forthcoming 
paper. 

4 RELIABILITY ANALYSIS OF AN 8-HOLE 
PLATE 

This structural reliability example is inspired from 
Deheeger and Lemaire ( |2007| ). It concerns the relia- 
bility analysis of a 200 x 100 mm 8-hole plate illus- 
trated in Figure |3j The diameter of the holes is set 
equal to 10 mm. Its left end is clamped both horizon- 
tally and vertically while its right end is subjected to 
a distributed line load with magnitude q = 100 MPa. 
Plain stress is assumed and the material is supposed 
to have a linear elastic behavior. The Poisson coeffi- 
cient 1/ is set equal to 0.3. Due to the boundary con- 
ditions the Poisson effect is not the same on all the 
plate though. The Young's modulus is modeled by an 
homogeneous lognormal random field with a mean 
^lE^ 200000 MPa, a coefficient of variation Se = 25% 
and assuming an isotropic squared exponential auto- 
correlation function with a 20 mm correlation length 
£. The two-dimensional random field is represented by 
a translated Karhunen-Loeve expansion discretized 
by means of a wavelet-Galerkin strategy proposed by 
Phoon et al.| P002[ ). The stochastic model involves 
20 independent standard Gaussian random variates 
grouped in the vector X to simulate the random field. 
The mechanical model is solved with Code_Aster 



( |eDF, R&D Division! , |2006| ) in order to retrieve the 
maximal Von Mises stress in the plate V. The perfor- 
mance function is then defined as follows: 



g{x) = (TO - niax{CTvon Mises (P)} 



(33) 



P/metalS ^PfeO^con 



(32) with respect to an arbitrary threshold ao = 450 MPa. 



(a) Mesh, boundary conditions and loads. 
Figure 3: Illustration of the 8 -hole plate example from 



(b) One sample path of the Young modulus random field. 
Deheeger and Lemaire|(|200'7]). 
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Table 1: Reliability analyses results for the 8-hole plate example from [Deheeger and Lemaire|p007| l. 



The proposed meta-model-based importance sam- 
pling procedure is applied to this structural reliability 
example. First an initial kriging predictor is built for 
the performance function g using a 100-point DOE. 
These 100 points are uniformally generated within 
the ^0 radius hypersphere. Based on this initial pre- 
diction, the DOE refinement procedure introduced in 



Section 2.3 is used. K = 100 new points are added at 
each refinement iteration. The refinement procedure 
is stopped after 1 000 estimations of the performance 
function. This may seem arbitrary but it is difficult to 
provide another stopping criterion for the refinement 
procedure - this needs further investigation. Then, the 
probabilistic classification function is defined with re- 
spect to the latest (finest) kriging prediction and it is 
used to compute the proposed estimator of the failure 
probability. 

The results are provided in Table [T] They are com- 
pared to a reference solution obtained by subset sim- 
ulation dAu and Beckl [200T]), and the multi-FORM 
estimator from Der Kiureghian and Dakessian ( 1998 1 
using FERUM v4.0 (Bou rinet et al.[|2009^ implemen- 
tations of these algorithms. FERUM is a Matlab tool- 
box for reliability analysis published under the Gen- 
eral Public License. The estimate of the augmented 
failure probability is equal to = 2.85 x 10~^, and the 
correction factor is equal to Scoit = 0.412. It means that 
the kriging predictor is rather accurate in that case. 
The probabilistic classification function is very close 
to its deterministic counterpart - and so is the instru- 
mental importance sampling density h. 

5 CONCLUSION 

Starting from the double premise that a surrogate- 
based reliability analyses does not permit to quantify 
the substitution error, and that the existing variance 



reduction techniques remain time-consuming when 
the performance function involves the output of an 
expensive-to-evaluate black box function, an hybrid 
strategy has been proposed. First, the probabilistic 
classification function was introduced, this function 
allows a smoother classification than its determinis- 
tic counterpart accounting for the epistemic uncer- 
tainty in the kriging prediction. Using this smoother 
classification function within an importance sampling 
framework then allowed to derive a meta-model- 
based importance sampling estimator. This estimator 
converges towards the theoretically impractical opti- 
mal importance sampling estimator and may provide 
a significant reduction of the estimation variance as 
illustrated in the example. 

In the present paper, the refinement procedure that 
leads to the probabilistic classification function is 
stopped arbitrarily. Work is in progress in order to es- 
tablish the best trade-off between the size of the DOE 
and the number of simulations required to estimate 
the correction factor acorr- 
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